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I.   INTRODUCTION 


The  parabolic  partial  differential  equation 


||  =  V (a(x,t,u) Vu)  +  f(x,t,u,Vu)     (x,t)  e  ftx(0,t] 


has  been  the  subject  of  the  development  of  several  numerical 
approximation  techniques.   The  equation  is  known  as  the 
heat  equation  in  a  specialized  form  and  is  useful  in  the 
study  of  heat  transfer  and  thermodynamics  in  general.   It 
also  has  applications  in  other  fields  of  engineering  and 
science,  such  as  fluid  dynamics  and  meteorology. 

The  development  of  Galerkin  methods  for  the  solution 
of  a  general  parabolic  and  hyperbolic  equation  has  been 
relatively  recent  [3,6,9].   This  thesis  describes  a  program 
for  the  solution  of  the  general  parabolic  problem  using  the 
central  different  Laplace  modified  alternating  direction 
Galerkin  scheme  described  in  [5,6,8].   The  program  was  built 
in  steps  beginning  with  the  one  dimensional  case  and  then 
extended  to  the  two  dimensional  case. 

The  program  solves  the  problem  for  a  rectangular  domain 
ft,  and  uses  nine  Gauss  quadrature  points  for  integration  in 
the  plane.   It  also  allows  for  nonlinearities  in  the  forcing 
function  f(x,t,u,Vu)  and  the  function  a(x,t,u)  can  be 
extended  to  have  nonlinear  terms  of  Vu.   The  theoretical 
proofs  of  stability  and  convergence  are  available  for  these 
cases  in  [4] . 


The  program  was  written  in  the  Fortran  IV  computer 
language  and  tested  on  an  IBM  360/67  computer  system. 


I I .   DEFINITIONS,  NOTATION,  AND  PRELIMINARY  CONCEPTS 

In  the  following  section  several  concepts  and  definitions 
which  are  essential  for  the  understanding  of  the  thesis  are 
presented.   Two  specific  spaces  of  functions  will  be  consid- 
ered and  definitions  of  the  norms  in  the  spaces  will  be 
given. 

Let  ft  c  R  ,  with  smooth  boundary  3°,,  be  the  domain  of 

interest  for  the  function  spaces.   The  first  of  these  spaces 

2 

is    that   of    square    mtegrable   or  measurable    functions,    L    (ft). 

2 

If      w  e L    (ft)  ,   then    the  norm  of  w   is   defined   as 


ILwHL2(pJ    =    £      w2    dx    <    °° 

o 

L  (ft)  is  a  Hilbert  space  with  respect  to  the  inner  product 


<f,g>  =  f   f(x)  g(x)  dx 

ft 

0  2 

The  Sobolev  space  H  (ft)   is  equivalent  to   L  (ft) . 

The  second  space  considered  will  be  the  Sobolev  space 
H  (ft) ,  defined  as 


H1  =  {w|DPw  £ L2(ft) ,  |p|  =  0,1} 


where 


D 


P_        8P 


|P|  =  %       (Pi) 


d   pl        9  Pn 

xl    * '  *    xn 
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The  norm  of  w  is  defined  as 


3w  .,  2     ,1/2 


HwIIh1(^)  =  [  'I W 1 1  X.2  C^)  +  ±l±   H33E7llL2(n,] 


Of  special  interest  for  the  development  of  this  paper 
is  the  space  Hn  a  subspace  of  H   defined  as 


H*  =  {w|w  eH1(fi)   and  w  =  0,  for  all  x  e  9ft} 


The  norm  for  this  subspace  is  defined  to  be 


ii  ,|         ,   _  .,  3w  |,  2     , 1/2 


A  major  concept  which  is  used  in  the  formulation  of 
Galerkin's  method  is  that  of  the  weak  solution  to  a  differ- 
ential equation.   Let  L  be  the  differential  operator 
describing  the  differential  equation,  that  is 

Lu  =  f  (2.1) 

A  function  u  is  said  to  satisfy  the  weak  form  of  equation 
(2.1),  if  for  every  test  function  v  in  some  test  space 

<Lu,v>  =  <f,v>  (2.2) 

If  u  is  not  sufficiently  dif f erentiable  to  be  substituted 
in  (2.1),  but  does  satisfy  (2.2),  then  u  is  said  to  be  a 
weak  solution. 
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To    illustrate    the   weak    solution,    consider 


Lu  =   u.   -  u        =    f(x,t)       ,  xel   =    (0,1),      t   >    0. 


u(0,t)    =   u(l,t)    =    0       ;       u(x,0)    =   uQ(x)     ,         xel 

(2.3) 


The  weak  form  of  this  equation  is 


<u.  -  u   ,v>  =  <f,v>     for  all  v£hJ(I)       (2.4) 


u(0,t)  =  u(l,t)  =  0 


u(x,0)  =  uQ (x) 


Integrating  by  parts  results  in 


<u.,v>  +  <u  ,v  >  =  <f,v>    for  all  veH:(I) 


u(x,0)  =  uQ(x)  (2.5) 


Note  that   ueL(I)   can  be  a  solution  to  (2.5),  that 

is  a  weak  solution  of  (2.3),  however  a  strong  solution  of 

2 
(2.3)  must  satisfy   ueC  (I)   for  each  t. 
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III.   THE  FORMULATION  OF  THE  GALERKIN  METHOD 

The  formulation  of  Galerkin  methods  for  the  solution  of 
differential  equations  is  based  on  the  weak  form  of  the 
equation.   Consider  the  nonlinear  parabolic  equation 


|H  =   V- (a(x,t,u) Vu)    +    f(x,t,u,Vu)       ,    (x,t)   eftx(0,T] 
u(x,0)    =   uQ(x)       ,      x  e  ft  (3.1) 

u(x,t)    =    0  (x,t)   e9ftx(0,T] 


and   the   test  space   HL(ft).      Clearly   for   any     veHQ(ft) 


•ttt  v  =   V* (a(x, t,u) Vu) v  +    f(x,t,u,Vu)v 

at 


Taking  the  integral  over  the  region  ft  yields 


/  |£  v  dx  =  /   V- (a(x,t,u) Vu) v  dx  +  /  f(x,t,u,Vu)v  dx 
ft  dr         ft  ft 


Using  Green's  first  formula  for  integration  by  parts  in 
multiple  dimensions  results  in 


/  -jnr  v  dx  =    /      a(x,t,u)Vu  v  dx  -    /      a(x,  t,u)  Vu*  Vv  dx 
ft    dt  ft  ft 

+    /      f(x,t,u,Vu)v   dx 
ft 
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Since  v  =  0  on  8ft  it  follows  that 


/   a(x,t,u)Vu  v  dx  =  0 

8ft 


and 


/  l£  v  dx  +  /  a(x,t,u)  Vu-  Vv  dx  =  /  f(x,t,u,Vu)v  dx 
ft  ft  ft 


Therefore  u  satisfies 


<y£  ,v>  +  <a(x,t,u) Vu,Vv>  =  <f (x,t,u,Vu) ,v>       (3.2) 
<u(x,0) ,v>  =  <uQ,v>       for  all  veHj(ft)  . 


i 

Let  M  be  a  finite  dimensional  subspace  of  Hfi(ft).   The 

continuous-time  Galerkin  method  is  to  find  for  each  t  e  (0,T] 
a  dif ferentiable  map   U(*,t):[0,T]  -»■  M  ,   such  that 


<~  rV>  +  <a(x,t,U)  VU,VV>  =  <f  (x,t,U,VU)  ,V>       (3.3) 

ot 


<U,V>  =  <uQ,V>  ,   t  =  0  ,   VeM 


Now  let  M  be  a  subspace  of  Hfi(ft)  such  that 

M  =  Span (w,,w2 / . . . ,w  )  ,   where   {w-j}-_i   1S  a  linearly 

independent  set.   Under  these  conditions  it  is  possible  to 

write 

n 
U(x,t)  =  l      a. (t)  w. (x)   . 
j=l  J  J 
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Substituting  in  equation  (3.3)  and  letting  V  =  w. 


n  n 

<  Z   a . (t)w. (x) ,w. (x) >  +  <a(x,t,U)   Z  a . (t) Vw. (x) , Vw. (x) > 

j=l  3  3  j=l  :     : 


=  <f (x,t,U,VU) /wi(x)> 


and 


n 


<  Z   a.(0)w.(x)>  =  <u0,wi(x)> 
j=l  •*  J 


Since  integration  is  a  linear  operation  and  the  inner  product 
defined  is  not  over  the  time  space,  it  follows  that, 


n    T  n 

Z   a^(t)<w. (x) ,w. (x)>  +   Z   a. (t) <a(x,t,U) Vw. (x) ,Vw. (x) > 

j=l   3      J  j=l  3  J 


=  <f (x,t,U, VU) ,w. (x) >  (3.4a) 


and 


n 


Z   a. (0) <w. (x) ,w. (x) >  =  <un/w.(x)>    i=l,2,...,n   (3.4b) 
j=l   3      D      i         u   i 


Equations  (3.4a)  and  (3.4b)  can  be  written  as  a  system 
of  n-nonlinear  ordinary  differential  equations 

Ma' (t)  +  S(a) a(t)  =  F 

(3.5) 
Ma(0)  =  b 
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where  S(a)  =  (s. .)  and  M  =  (m. .)  are  n    order  matrices  and 


s. .  =  <a(x,t,U) Vw. (x) , Vw. (x) > 
13  J       1 


mii  =  <wi  ^  /Wi  ^  > 


and 


F  =  (fi)  ;    f±   =  <f (x,t,U,VU) ,wi(x) > 


T 
a  =  (a,  ,a2 / . . . r a  ) 


T 
b=  (<uQ/w1>/<u0,w2>/ . . . /<uQ/wn>) 


Using  the  properties  of  the  functions  w. f  it  can  be 
shown  that  the  matrix  M  is  positive  definite,  and  with  the 
additional  property  that  a(x/t/u)  is  bounded  above  and  below, 
the  matrix  S  can  also  be  shown  to  be  positive  definite.   If 
in  addition  a(x,t,u)  and  f(x,t,u,Vu)  satisfy  Lipschitz 
conditions,  it  can  be  shown  that  the  system  of  differential 
equations  has  a  unique  solution. 

In  the  one  dimensional  case,  an  approximate  solution  to 
equation  (3.5a)  can  be  obtained  by  combining  it  with  a 
finite  difference  scheme  to  discretize  in  time  as  follows. 
Assuming,  for  simplicity,  that  a(x,t,u)  =  a  =  constant 
and  f(x,t,u,ux)  =  f(x,t) 

m+1   m       m+1  ,  m       ,  , 
M(a   At~a  )  +  S(a   2+a  )  =  Fm+*  (3.6) 
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This  is  a  system  of  n  equations  and  n  unknowns  and  can  be 
simplified  to 


(2M+AtS)am+1  =  (2M-AtS)am  +  2AtFm+^ 


or 


(2M+AtS)  (am+1-am)  =  -2AtSam  +  2AtFra+^ 


where   am  =  a.  (t  )   and  F™+k  =    <f(t  .J,w.> 
j     3   m         1         m+%  '  1 


In  the  above  equation  the  time  variable  t  has  been  discretized 
to   t  =  mAt  ,   where   t  =  T/N  ,   N  a  positive  integer. 
This  method  is  known  as  the  Crank-Nicholson-Galerkin  method. 

The  one  dimensional  program  was  written  to  set  up  the 
matrices  defined  above  and  solve  the  algebraic  equation  (3.7). 
A  minor  alteration  in  the  program  would  allow  the  use  of 
(3.7a),  and  give  the  advantage  of  reducing  the  computational 
roundoff  error.   The  program  is  also  limited  to  the  linearity 
simplifications  made  above. 

The  required  accuracy  of  the  method  in  the  two  dimensional 
case  leads  to  a  large  number  of  basis  functions  for  the  space 
and  consequently  would  require  excessively  large  matrices. 
The  matrices  have  a  band  structure  and  the  band  width  of  the 
matrices  is  dependent  on  the  number  of  intervals  used  in  the 
grid  set  up  for  the  region.   These  large  matrices  can  be 
avoided  by  implementing  an  alternating  direction  method,  and 
using  a  tensor  product  basis  for  the  space.   This  tensor 
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product  basis  reduces  the  problem  to  working  with  matrices 
which  have  a  constant  band  width  of  seven,  in  the  case  of 

M  being  Herrnite  cubic  splines.   It  also  greatly  simplifies 

2 
the  problem  when  the  region  U     R   is  rectangular. 

Let   M   =  Span (£, (x) , £~ (x) , . . . , £   (x) )   be  a  basis  for 
x         i      z         n 

x 

one  dimension  of  the  space  M  and   M   =  Span(n, (y) f n~ (y) , . . , 

y       i    z 

n   (y) )   be  a  basis  for  the  other  dimensions  of  M.   M   and 

ny  , 

M   are  bases  for  one  dimensional  subspaces  of  Hrt(I  )  and 
y  r  0   x 

Hq(I  ),  such  that   M  =  Mx  fi  M   =  Span  (  ^r^  ,  £  n2  '  '  '  '  '  ^lnn  ' 

y 

. . . , £   n   )  .   Also  define 
'  n   n 
x  y 


<f,g>   =  /   fg  dx 

X    I 
x 


<f,g>,r    /   fg  dy 

y 


y    i 


Since   M  SM    forms  a  basis  for  the  space  M,  we  can  write 
x    y  tr  i 

n   n 
x   y 

U(x,y,t)  =11        a   (t)(L(x)8H  (y))  . 
p=l  q=l    P4      P       4 


There  are  several  time  discrete  methods  with  which  to 
proceed  at  this  point.   The  program  has  been  designed  to  use 
the  alternating  direction  version  of  the  centered  difference 
Laplace  modified  Galerkin  method  described  in  [5,6,8]  and 
defined  as 
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m+1  _  TTm-l 
<  - — 2At^ /V>  +    <a    te't'W  VUm,  W> 

+  2x2At<"i  ^  4^>  =  <fm<*,t,u,vu),v> 


and 


<U,V>  =  <uQ,V>     VeMxfiMy  (3.7) 


where  the  superscript  of  the  functions  denotes  the  time 
step  at  which  it  is  to  be  evaluated,  and 


92Um  =  Um+1  -  2Um  +  Um_1 


and 


32V    T2,m 
£  L  (ft) 


dxdy 
X   is  subject  to  the  restriction 

X    >  4  amax  =  4"  maX  a  ^X' t  ,u^x' t)  * 

x  ft 

t>0 

for  stability  [5,6]. 

Equation  (3.7)  can  be  written  in  the  algebraic  tensor 
product  form 

(C  +  2AAtA  )  a  (C  +  2AAtA  )  am+1  =  2(2AAt(C  ®A  +A  fiC  ) 
x       x     y       y  x   y   y   y 

+  4X2(At)2(AvSA  ))am 

+     (CxQCy-2AAt(CxSAy+AxQCy)   -  4  A2  (  At)  2  (Ax  fi  A    )  )  am_1 

+    2At$m 
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where 


«pq  =  -<am(x,t,u)vum,vup  nq)>  +  <fm(x,t,u,vu) ,£p  nq> 


and 


x    *i   j  x  y     i   j  y 


A   =  <£.,£.>  A   =  <n?",n*> 

x    %i'*3  x  y     'i'  'j  y 


t-n  e  ■      ni    m   m-1     .,        .  .    ,_  _.      , 
Define   e   =  a  -a     ,   then  equation  (3.7)  can  be 


written  as 


(C  +2AAtAj  Q  (C  +2AAtA  )  (em+1-em)  =  -2(C„fiC  )em  +  2At$m 
x       x      y       y  y 


and 


(3.8a) 


m+1    ,  m+1   im  ,   m  .   m  , ,,  n,  x 

a    =  (e   -e  )  +  e   +  a  (3.8b) 


The  scheme  (3.8)  adds  significant  efficiency  to  the  program 
and  also  reduces  roundoff  error. 

The  Laplace  centered  difference  method  initially  requires 
two  sets  of  the  coefficients,  those  at  time  zero  and  at  one 

time  step  later.   At  time  zero  the  set  of  coefficients  are 

2 
obtained  by  an  L  (Q)    projection  of  the  initial  conditions 

into  the  space,  but  the  second  set  presents  a  problem  where 

a  variety  of  methods  are  available.   Two  of  these  methods 
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were  attempted,  expanding  the  initial  conditions  in  a 
Taylor's  Series  about  t  =  0  and  using  the  Laplace  forward 
difference  method.   The  former  involved  taking  the  deriva- 
tives of  the  initial  condition  function  and  presented  a 
problem  when  these  derivatives  resulted  in  zero.   The  method 
which  was  finally  used  was  the  Laplace  forward  difference 
scheme.   It  is  defined  as  follows 


(C  +  XAtA  )  8  (C  +  XAtA)  (am+1  -  am)  =  (At)  4>m      (3.9) 
a       rft      y        y 


At  first  glance  the  matrices  in  (3.8a)  and  (3.9)  seem  to  be 
different,  however  in  the  case  of  this  method  X  is  under 
the  restriction 


X  >  yr   a    =  yr   max  a(x,t,u)     x  e  0.    ,      t>0 
2   max    2 


thus  a  prudent  choice  of  X  results  in  identical  matrices  and 
it  is  not  necessary  to  form  a  second  set  of  coefficient 
matrices  in  the  program.   It  is  also  noteworthy  that  <P   has 
the  same  form  as  in  the  previous  method  and  the  subroutine 
used  to  compute  $  may  be  used  for  both  with  a  correction  of 
y  for  the  forward  difference  scheme.   These  advantages  led 
to  the  final  choice. 

The  tensor  product  of  matrices  has  the  property  that 


(B  QD)  =  (BSD  (I  SD) 
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This  property  enables  a  significant  simplification  of  two 
aspects  of  the  programming  problem.   First,  the  tensor 
product  multiplication  on  the  right  hand  side  of  the  equation 
can  be  accomplished  by  computing  first 


.  _  _  ^  x  m    m 
(I. fit  C  )e  =  y 


and  then 


(CxSI)ym 


By  properly  ordering  the  vector  em ,  this  results  in  a  series 

of  multiplication  by  the  relatively  small  matrices  C   and  C  , 

x      y 

that  is,  matrices  which  have  a  constant  band  width  of  seven 
for  Hermite  cubics  rather  than  a  band  width  which  depends 
on  the  number  of  intervals  in  the  grid  of  the  problem  as  in 

c  a  c  . 

x   y 

This  same  idea  can  be  used  in  solving  the  linear  system. 

Defining  the  right  hand  side  of  the  equations  (3.8,3.9) 
to  be  $  ,  let 


(I(MCy+2AAtAy))(em+1-em)  =  Y™ 


Using  this,  solve  first 


((C     +2AAtA    )  Sl)ym  =    3m 


and   then    solve 


(IS  (Cy  +  2XAtAy))  (em+1-em)    =ym 
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Again  proper  ordering  of  the  vectors  involved  results  in 
the  simplification  of  solving  a  system  of  one  dimensional 
problems  rather  than  the  larger  problem  of  solving  the  full 
tensor  product  matrix. 

The  two  dimensional  program  has  been  designed  to  solve 
the  problem  in  the  manner  described  above.   The  basis 
functions  used  in  the  program  are  Hermite  bicubic  splines, 
taken  as  the  tensor  products  of  one  dimensional  Hermite 
cubic  splines,  and  will  be  described  in  Section  IV.   The 
procedure  described  above  has  the  drawback  that  it  is  in 
general  good  only  on  rectangular  polygons,  and  must  be 
altered  significantly  to  handle  regions  other  than  simple 
rectangles.   It  can  be  changed  to  handle  regions  which  are 
unions  of  rectangles  [5] . 
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IV.   HERMITE  BICUBIC  SPLINES 

The  Galerkin  procedure  is  dependent  on  a  set  of  linearly 
independent  basis  functions  for  the  space  considered  in  the 
problem.   Hermite  bicubic  splines  provide  such  a  basis. 
These  splines  are  a  high  powered  tool  for  the  approximation 
of  smooth  functions. 

The  Hermite  bicubic  splines  are  taken  as  the  tensor 
product  of  one  dimensional  Hermite  cubics.   These  cubics 
have  a  piecewise  polynomial  nature,  where  the  polynomial  is 
of  degree  less  than  or  equal  to  three.   Hermite  cubics  have 
globally  continuous  first  derivatives.   Thus,  the  tensor 
product  of  the  one  dimensional  cubics  has  the  property  that 
both  first  partials  and  the  first  mixed  partial  derivatives 
are  continuous.   These  properties  are  sufficient  for  the 
approximations  desired  as  a  result  of  the  program  [1] . 

The  package  of  subroutines  described  in  [2]  was  used  for 
the  program.   It  consists  of  subroutines  which  compute  the 
values  of  the  B-splines  at  specific  points,  and  given  a  set 
of  coefficients,  computes  the  value  of  the  approximated 
function  at  specified  points.   The  high  degree  of  flexibility 
of  the  package  was  the  major  reason  for  its  use  in  the 
program.   The  program  can  be  extended  to  use  higher  degree 
B-splines  which  gives  a  greater  degree  of  smoothness  and 
accuracy  to  the  approximating  function.   Also,  the  extension 
of  the  program  to  higher  dimensions  can  be  accomplished  by 
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taking  a  higher  degree  tensor  product,  and  making  the 

appropriate  alterations  in  the  computational  procedures. 

Consider  the  interval  [0,1],  and  the  partition  of  the 

interval  A:   0  =  x^  <  x,  <  ...  <  x  =  1  .   Let   h.  =  x.  -x.  ,, 

0—1—     —  n  11    l-l 

and   I.  =  [x.  ,,x.].   Define   h  =  max  h.  ,   and   h  =  min  h.  . 
l      l-l'  i  .    i         —     .    l 

The  Hermite  cubics  are  piecewise  cubics  over  I.  for  all  i, 
and  are  members  of  C  (I) .   The  order  of  accuracy  for  Hermite 
cubics  is  given  by  the  following  well-known  theorem: 

4 
Theorem.   If   u  e  H  (I)   and   h/h  <  a  <  °°  ,   then 


inf  ||u-u||  hS(i)  <  Ch4  S  ||u||  h4(i)   ,   0  <  s  <  3 


where  M  is  the  space  of  Hermite  cubics. 


25 


V.   ERROR  ESTIMATES 

The  usefulness  of  a  procedure  to  obtain  an  approximate 
solution  to  a  differential  equation  is  measured  by  its 
degree  of  accuracy  and  the  efficiency  and  costs  of  its 
implementation.   The  purpose  of  this  section  is  to  derive 
an  error  bound  for  the  continuous  time  Galerkin  method  for 
a  linear  version  of  equation  (3.1).   The  following  theorem 
and  its  proof  is  found  in  [8],   It  is  restricted  to  the  case 
where  a(x,t)  is  not  a  function  of  u.   Proofs  for  the  nonlinear 
case  are  also  available  [3,6,9]. 

Theorem  1.   Let  u  be  the  solution  to 


<lr  'v>  +  <a(x't)  Vu,Vv>  =  0    vehJ(^)  ,   te  (0,T] 


(5.1) 


<u(x,0)>=  <uQ(x)>    x  e  ft 


and  let  U  be  the  solution  to 


<-j~    ,V>  +  <a(x,t)  VU,  VV>    =    0  VeM    ,       te(0,T]       (5.2) 


<U,V>    =    <uQ,V>  t   =    0 

and 

0    <    c0    <    a         (x,t)     <    c,  (5.3) 

0   —     max  —      1 
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Then  there  exists  constants  C  and  9,  which  depend  on 
cn,  c, ,  and  T,  such  that 


KMT   HU-<2(n)(t)  +  3  /  l|u-U||2,(t)  dt 


<  c(0^T  ||u-a||2a(a)(t)  +  /  l|u-a||gi(0)(t)  at 

+  /  ||  |^  (u-u)||22(n)(t)  dt)        (5.4) 


where  u  is  an  arbitrary  map  of  (0,T]  into  M. 

Proof:   Replace  v  with  V  in  (5.1)  and  subtract  (5.2)  from 
(5.1) 


<  ~    ,V>  +  <a(x,t)Vu/VV>  =  0 
dt: 


<  %rr   rV>   +  <a(x,t) VU,VV>  =  0 

a  u 


results  in 


<  ~t  U)  'v>  +  <a(x,t)  V(u-U)  ,VV>  =  0   . 

dt 


Let  e  =  u-U   and  let  V  =  e+e  ,   where   e  =  u-u  ,   then 


<  |§  ,e+e>  +  <a(x,t) Ve, Ve+e>  =  0 

o  t 
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Note 


<  at  'e>  =  £    Ft  e  dx 


1   d         r    ~2     , 
^-  ^rr      /    e      dx 

2  dt    n 


i^-     Hell  29 
2    dt      l|e|1  L2(ft) 


and 


and 


<aVe/Ve>    >    cQ   ||e||  21(fl)  from    (5.3) 


|<aVe,Ve>|    <      ||aVe||  L2(fi)     -      ||Ve||  L2  ( 


ft) 


1  Cl  HellHj(ft)  *   H^H^) 


from  Schwarz ' s  inequality. 

Young's  inequality  states   ||f||  •  ||g||  <  e  ||f|    +  —  ||g 
Using  this  inequality  we  get 


<ave,ve>|  <  eC;L  ||e||  2i(fi)  +  47  ||S  ||  *!  (pj 


Thus 


li.      HpII  2o  +    <le-      e>    +    c         llell  2i 

2    dt      llell  L2(ft)  <8t    'e>    +    C0      llell  hJ(«) 


2  Cl  2 

<    cc,      llell  „i/n.    +  —     llell  „i 


Hj(ft)  4c      llc"  Hj(fi)       t  e  (0,T] 
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and 


Id        i,     i,  2 
ii  e  ii 


8e 


>L2(fi)     +    <It    ,e>    +    (C0      ecl)    lle'lHl(fi) 


1         ~      2 

-  47    Hell  hJ(^)      ' 


t  e  (0,T] 


Let      e    =    0^/20,     ,       then    for      C   =    c,    /c 
0         1  1        o 


"eH  L2(fi)    +    2    <!!    >*>    +    c0   il°H  Hl(P0 


dt 


^  c  "-ii  Hi(Q) 


Integrating  from  0  to  t 


e"£2(n)<T)    '      'HIP 


{fl)(0)    +   cQ    /        ||e||Hl(n)    dt 


and 


±c  0/T  hsHhJ(«)  dt  - 2  ^  41  'g>  dt 


/      <!f   ,e>   dt  =   <e,e> 
0         dt 


-  {T  <e'  H> dt 


<e,e>|    <     e||e||  22(n)    +   C  ||S!l  ^(Q) 
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and 


/<e,ff|    dt  <     \   (/  ||e||  22(fi)  (t)    dt  +   /   ||||||  22(fi)  (t)    dt) 


It  follows  that 


2  2  T  0 

0  H0 


e'l  L2(fi)  CT)     "      lleIU2(fi)(0)     +    C0    £    IMlHitM     dt 


<     ||e||22(n)(T)    +C(   ||e||22(n)(0)    +   /T||e|l^l(n)    ^ 

+   flMI&Cfl)    dt+      H§llL2(n)(T)    +      HgllL2(fi)(0) 
+  0/TH|fllL2(«)    dt) 


Thus  for  e  sufficiently  small 


HellL2(fi)(T)  +  3  0'THellHl(«)  dt 

1    C  (0^llellL2(n)  dt+   HellL2(n)(°»  +  jf  PHh^)  dt 
+   ||e||^2(n)(T)  +  l|g||^2(n)(0)  +  /||f||22(fi)  dt) 


Gronwall's  inequality:   Let  f,g,  and  h  be  piecewise  continuous 
nonnegative  functions  defined  on  an  interval   a  <  t  <  b   and 
assume  that  for  each  t  e  [a,b] 


T 

f(t)  +  h(t)  <  g(t)  +  /   f(s)  d! 

0 


30 


Then 


r   t-s 
f (t)  +  h(t)  <   /   e"   g(s)  ds   +   g(t) 


and  if  g  is  nondecreasing 


f(t)  +  h(t)  <  efc  a  g(t)  . 


Thus 


•\\hla)M  +  *  (    IIHl£i{0,   dt 


iC    <0<?<T      llellL2(fi)(t)     +    /    IIS  11^1(0,     dt 


+    (    liffH  L2(„,     «+      HSHlW0)     > 


Now  by    the    definition   of    U(x,0) 


u-ullL2(fi)«»    1     l|u-0||    L2(Q)(0) 


and  the  desired  result  is  obtained,  that  is 


o<t<T    \\^\\hUl)^  +  a  /    ll*-n||2J(fl)   at 


C    <0<?<T      ll«-«ll£»(0)<*>     +    ^    llU-aHHl(fi)(t)     dt 

+ /lift  (u-a)||22(n)(t)  at  ) 
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This  implies  that  the  error  bounds  are  dependent  on  how  well 
u  and  its  first  partial  derivative  with  respect  to  time  can 
be  approximated  by  the  arbitrary  map  u.   For  Hermite  cubics, 
||u-u||L2(fi)   is   0(h4)  ,   and  l|u-u||Hl(n)   is  0(h3)  . 

With  further  restrictions  it  can  also  be  shown  that  the 

2  4 

L  (fi)  error  is  0(h  ) . 

In  [3]  J.  E.  Dendy,  Jr.  proves  the  following  theorem 

which  gives  error  estimates  for  the  discrete  time  case. 

Theorem  2 .   Let  u  be  the  solution  of 


<ft  '    v>  +  <a( ' /U)  Vu'Vv>  =  <f(*/U),v>  ,    veHjJ(fl) 


u(x,0)  =  Uq (x)       x  e  tt 


u(-,t)   Hq(^)1   ,    t  e  (0,T) 


and  assume  that   A  >  j   c,  ,  c1 =  max  a(x,t,u) , (x,t)efix(0,t] . 

Let  Em  =  u^-W™  ,   where  Um  e  M,   and  M  satisfies 

inf  ||u-u||H3(Q)  1  C  hP"*S  ||u||Hr(Q)  [3].   Let  W  be  the 

mapping  from:   (0,T]  +  M  defined  by 

<a(-,u  (u-W)  ,V>  +  <u-W/V>  =  0   for  all  VeM  . 

2 

If   ueC4(fiM0,T]),  and   ^Ip^T  ^ (0  'T;Hr  (fi)  )  '   r-2' 

a  t 

then  for   1  <  N  <  M   and   h  and  At  sufficiently  small 
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N-l 

~       ~      ii^_,mn  2     .,  ,    |i_Nii  2  n„N-lii  2, 

8       E      II  9E     ||        At   +    q[      E      |   .    +         E  ||   n  ] 

m=l  L  X 

+    2A2At    ||  |V    (EN-EN-1)||2    <    eCAtN(At)4 
11   dxdy  '  "        — 


N    MN.  ,  |  cNAt/A4_.2  cNAt   .  r 

u  -U    )  |       =   e  (At)       =   c   e  h         , 


if 

HE"!     J    +      II E"    |  l   +    At    |-^    (EA-E0)  ||  -    <    c(At) 


/T,1      ^,1,2^       /A^N  4 


For  the  case  of  Hermite  cubics ,  r  =  4,  and  the  conditions 
on  M  are  satisfied.   U  is  the  approximate  solution  obtained 
in  the  Laplace  modified  alternating  direction  method. 

Using  the  program  it  can  be  shown  that  the  order  of  the 
error  is  actually  this  by  approximating  the  value  of  oo  where 
the  error  is  CKh^).   This  approximation  can  be  obtained  in 
the  following  way:  " 

Let  e,  be  the  error  for  h  =  h, ,  and  e~  be  the  error  for 
h  =  h2  r      then 


e   =  C  h£  e2  =  C  h^ 


log(e1/e2)  =  oo  log(h1/h2) 


0) 


=  (log(e1)  -  log(e2))/(log(h1)  -log(h2)) 
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1/2 

Using  this  method  and  choosing   t  =  h     we  obtain  the 

results  presented  in  Table  1.   Table  2  shows  similar  results 
for  the  time  error. 

To  obtain  the  estimates  in  the  two  tables  the  following 
problem  was  used  for  easily  obtained  analytic  solution. 


ut  =  uxx  +  uyy  -  2(e  fc) ((x-x2)  +  (y-y2)  +  |(x-x2) (y-y2) ) 


u(«,t)  =  0      for  all  x£3fi 


2      2 

u(x,y,0)  =  (x-x  ) (y-y  ) 


In  this  test   A  =  -r     was  used  in  the  discrete  time 

4 

formulation . 


TABLE  1 
H  DT  L-2  NORM        MAX  NORM 


0.2000      0.0400         0.0003680       0.0006815 

0.1000      0.0100         0.0000286       0.0000619 

THE  VALUE  OF  W  FOR  THE  MAX  NORM  IS  3.46006. 
THE  VALUE  OF  W  FOR  THE  L-2  NORM  IS  3.68510. 

TABLE  2 
H  DT  L-2  NORM        MAX  NORM 


0.2000      0.0400         0.0003680       0.0006815 

0.1000      0.0100         0.0000286       0.0000619 

THE  VALUE  OF  W  FOR  THE  MAX  NORM  IS  1.7  300  3. 

THE  VALUE  OF  W  FOR  THE  L-2  NORM  IS  1.84255. 
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VI.   TEST  PROBLEMS 

The  utility  of  the  program  is  evident  in  several  fields 
of  engineering.   The  following  test  problems  will  illustrate 
its  use  in  two  of  these.   The  problems  are  taken  from  the 
fields  of  fluid  dynamics  and  heat  transfer. 

A.   RAYLEIGH'S  PROBLEM  FOR  A  CORNER 

In  the  first  problem  we  consider  the  diffusion  equation 
in  two  dimensions  with  the  following  boundary  and  initial 
conditions 

ut  =  (1/Re)  (uxx  +  Uyy)     ,    x,y  >  0 

u(x,y,0)  =  1 

u(0,y,t)  =  u(x,0,t)  =  0 

u(x,y,t)  ■*•  1   as   x  and  y  ■+  °° 

where  Re  is  the  Reynold's  number  of  the  fluid  under 
consideration. 

This  is  known  as  Rayleigh's  problem  for  a  corner,  and 
the  solution  describes  the  impulsive  motion  of  a  right 
angled  corner  formed  by  two  infinite  flat  plates  and  is 
used  to  infer  the  steady  flow  along  the  corner,  with  leading 
edge  at  t=0 . 


35 


For  the  purpose  of  numerical  approximations,  using  a 
Reynold's  number  of  1000,  x  =  4  and  y  =  4  can  be  considered 
as  infinity  and  the  problem  can  be  stated  for  use  in  the 
program  as 


u.  =  (1/1000)  (u   +  u   ) 
t  xx    yy 


and 


u(x,y,0)  =  1 

u(0,y,t)  =  u(x,0,t)  =  0 

u(4,y,t)  =  u(x,4,t)  =  1 

a(x,y,t,u)  =  1/1000 


f(x,y,t,u,ux,u  )  =  0 


A  nonuniform  grid  with  a  fine  mesh  near  the  x  =  0,  y  =  0 
corner  is  useful  for  an  accurate  description  of  the  boundary 
layer  behavior  near  the  corner.   Illustration  (1)  is  a 
diagram  of  the  grid  used. 

The  exact  solution  to  this  problem  can  be  found  analyt- 
ically to  be 

u(x,y,t)  =  erf(X)  erf(Y) 

where 

X  =  (x/2) (Re/t)1/2 

Y  =  (y/2) (Re/t)1/2 


36 


1 

1 

— 

— 

u 

I 

,    . 

1 

.  . 

Illustration    (1) 


The  nature  of  the  solution  is  that  of  a  shock  wave 
traveling  through  the  rectangle.   The  error  in  the  approxi- 
mation seems  to  be  greatest  on  the  wave  and  small  even  when 
the  grid  size  is  quite  large.   The  error  was  also  larger  at 
early  times  near  the  boundary  discontinuity  than  later  in 
the  process.   Table  3  is  a  sample  of  the  output  of  the 
program  for  this  problem.   The  output  points  were  taken  near 
the  discontinuity  on  the  boundary  since  it  is  this  area 
which  is  of  greatest  interest. 

For  this  problem   A  =  .00025,  At  =  .03  ,   and  the 
largest  Ax  was  .75. 

B.   THE  HEAT  TRANSFER  PROBLEM 

The  equation  for  heat  conductance  for  heterogeneous 
isotropic  solids  and  frictionless  incompressible  fluids  is 


pc  |H  =  V(k  u)  +  q0'" 


where  the  solution  u(x,y,t)  represents  the  temperature  of 
the  material.   p  is  the  density  of  the  material,  c  is  the 
heat  capacitance,  k  is  the  thermoconductivity ,  and  qQ1"  is 
an  internal  generation  term.   pc  is  a  constant  and  can  be 
combined  with  k  and  q0,M  to  give 


|H  =  v(k'  u)  +  q-'/pc 


In  this  form  k1  is  the  thermal  dif f usitivity 
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For  most  materials  the  thermal  dif fusitivity  can  be 
approximated  by  a  linear  function  of  the  temperature.   It 
is  typically  of  the  form 


k'  =  k  +  a(u-uQ) 


where  k  is  the  initial  value  of  the  thermal  dif fusitivity 
and  a  is  a  small  real  number. 

The  internal  generation  term  is  highly  variable  and 
dependent  on  the  type  of  problem.   A  typical  function 
encountered  is  of  the  form 

e$(u/u0) 

Analytic  solutions  for  problems  such  as  this  are  not 
available  and  it  is  in  these  cases  that  approximation 
techniques  are  invaluable. 

The  problem  for  the  program  run  will  be  as  follows 


ut  =  V  (  (100  +  .005(u-u0) ) Vu)  +  e(u/uo) 


u(x,y,0)  =  sin(iTx)  +  sin(Try) 
u(0,y,t)  =  u(l,y,t)  =  e{~v'    t)  sin(Try) 
u(x,0,t)  =  u(x,l,t)  =   eK    ^      '    sin(TTx) 


a(x,y,t,u)  =  100  +  .005(u-  (sin(iTx)  +sin(iTy))) 

N     (u/(sin(iTx)+sin(7Ty)  )  ) 
f (x,y,t,u,u  ,u)  -  ev 
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with  uniform  grid,  with   Ax  =  Ay  =  h  =  1/10  ,   time 
step  At  =  .01   and    X  =  30  . 

Plots  1-8  show  the  general  nature  of  the  solution. 
Plot  1  shows  the  solution  at  an  early  time  while  it  retains 
much  of  its  initial  form.   The  area  near  the  boundary  is 
beginning  to  become  smaller.   The  internal  generation  term 
causes  the  oscillatory  nature  of  the  solution.   After  a 
longer  period  of  time  the  solution  begins  to  take  on  a 
steady  state  of  zero. 
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PLOT    (1)        TIME    =    .05 
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K 


PLOT     (2)       TIME    -     .40 
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PLOT    (3)       TIME    =    .5 
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PLOT     (4)       TIME    =    .6 


PLOT    (5)       TIME    =    .7 
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PLOT     (6)       TIME    =    .9 
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PLOT    (7)       TIME    =     .11 
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PLOT    (0)       TIME    =    2.5  5 


VII.   CONCLUSIONS  AND  RECOMMENDATIONS 

Two  major  conclusions  can  be  reached  from  the  applica- 
tion of  the  program  to  the  various  problems.   First  the 
program  is  in  general  more  efficient  than  the  standard 
implementation  of  Galerkin  methods  which  achieve  the  same 
degree  of  accuracy.   It  is  also  possible  to  achieve  a  higher 
degree  of  continuity  and  approximation  in  the  solution  using 
higher  degree  B-splines.   It  is  possible  that  changes  not 
evident  to  the  author  may  be  made  to  improve  the  efficiency 
of  the  program.   Most  of  the  computational  effort  in  a  given 
time  step  is  in  evaluating  U   and  VUn  at  the  gauss  points; 
more  efficient  evaluations  or  placement  of  the  quadrature 
points  can  possibly  increase  the  speed  of  the  program.   It 
should  be  noted  that  the  matrices  [see  (3.8) ,(3.9)]  arising 
need  be  factored  only  once  since  they  are  time-independent. 
Second,  optimal  error  bounds  for  the  method  have  not  been 
derived  in  the  case  that  the  Laplace  modified  forward 
difference  method  is  used  to  obtain  the  second  set  of 
coefficients  so  that  the  centered  difference  method  may  be 
used,  but  this  does  not  seem  to  have  any  effect  on  the 
solution  in  application. 

Using  the  program  as  a  basis  the  Galerkin  methods  may 
be  extended  to  a  third  dimension,  using  extension  of  the 
alternating  direction  to  obtain  it.   The  basic  ideas  used 
in  the  development  of  the  program  can  be  utilized  to  program 
the  second  order  hyperbolic  problem  using  the  Galerkin  methods 
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APPENDIX  A 
THE  COMPUTER  PROGRAM 

A.  THE  PROGRAM  VARIABLES 

Alpham  -  the  coefficients  of  the  splines  at  time  m. 
CX,CY,AX,AY  -  the  matrices  C  ,C  ,A  ,A   respectively 

X   y   X   y 

described  in  Section  III 
X,Y  -  the  grid  points  on  the  x  and  y  axes 
TKX,TKY  -  the  knots  system  on  the  x  and  y  axes,  used 

for  the  de  Boor  routines 
GX,GY  -  the  Gauss  weights  computed  in  WEIGHT 
EKCX,ETAY  -  the  values  of  the  basis  functions  at  the 

quadrature  points 
ALXR,ALXL,ALYB,ALYT  -  the  coefficients  for  the  projection 

of  the  boundary  conditions  at  time  m 
CPX,CPY  -  the  quadrature  points  on  the  x  and  y  axes 
NX, NY  -  the  number  of  intervals  on  the  x  and  y  axes 
IAX,IAY  -  the  dimension  of  the  matrices  defined  above 
DT  -  At 
TIME  -  time  m 

B.  THE  MAIN  PROGRAM 

The  main  program  calls  various  subroutines  to  accomplish 
the  necessary  procedures.   Brief  indications  of  the  function 
of  each  subroutine  will  follow.   The  program  begins  by 
setting  up  the  grid  desired  and  its  corresponding  knots 
system  and  then  computes  the  necessary  matrices.   It  stores 
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parts  of  the  matrices  which  will  be  changed  but  needed  for 
future  corrections.   The  first  set  of  coefficients  is 
obtained  by  projecting  the  initial  conditions  into  the 
space.   The  second  set  is  computed  using  a  forward  differ- 
ence Laplace  modified  step.   The  program  then  steps  forward 
in  time  using  the  central  difference  Laplace  modified 
method. 

C.  MAT COM 

MATCOM  computes  the  matrices  for  the  program  and  stores 
them  in  band  symmetric  storage  used  by  IMSL  library  routines 
It  also  stores  the  quadrature  points  and  the  values  of  the 
basis  functions  at  the  quadrature  points. 

D.  FMATCO 

FMATCO  projects  the  initial  conditions  into  the  space. 

E.  WEIGHT 

WEIGHT  computes  and  stores  the  product  of  the  quadrature 
weights  and  Ax  and  Ay  for  the  appropriate  interval. 

F.  GRID 

GRID  computes  the  grid  points. 

G.  KNOTS 

KNOTS  computes  the  knot  system  for  the  de  Boor  routines. 

H .   UEVAL 

UEVAL  computes  the  value  of  the  function  at  specific 
points  (because  of  inefficiency  its  use  should  be  limited) . 
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I.   MULTIP 

MULTIP  computes  (CX   CY)  em   in  manner  described  in 
Section  III. 


J.   SOLVER 

SOLVER  solves  system  of  the  form  (CX   I) (I   CY)  em  =  3 
in  the  manner  described  in  Section  III. 


K.   ANEVAL 

ANEVAL  evaluates  the  function  and  its  first  partial 
derivatives  at  the  quadrature  points. 

L.   RTHDSD 

RTHDSD  computes  the  right  hand  side  of  equation  (3.8). 

M.   ALXYBC 

ALXYBC  projects  the  boundary  conditions  into  the  space, 

N.   FUNCTION  F 

Function  F  defines  the  necessary  functions  for  the 
problem  being  solved. 

0.   THE  DE  BOOR  SUBROUTINES 

The  de  Boor  subroutines  perform  various  operations  to 
obtain  B-splines  [2]. 
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APPENDIX  B 
THE  L2  PROJECTION  OF  FUNCTIONS 

2 

The  L   projection  of  the  initial  conditions  and  other 

functions  involved  in  the  program  were  accomplished  using 
Gauss  quadrature.   Appendix  A  describes  several  subroutines 
in  which  this  is  done.   To  achieve  the  desired  accuracy  it 
was  necessary  to  use  nine  quadrature  points  for  each  interval. 
In  the  subroutine  WEIGHT,  the  Gauss  weights  are  combined  with 
the  value  of  Ax  and  Ay  for  their  respective  intervals  and 
saved.   This  is  necessary  in  the  case  of  non-uniform  grids. 
The  quadrature  points  are  stored  for  future  use  in  the 
subroutine  MATCOM. 

To  project  a  function  into  the  space,  we  must  take  its 
inner  produce  with  each  of  the  basis  functions.   This  means 
we  must  know  the  values  of  the  basis  functions  at  the 
quadrature  points.   These  are  computed  using  the  de  Boor 
routines  and  saved  in  MATCOM.   In  the  case  of  one  dimensional 
B-splines  there  are  only  four  non-zero  basis  functions  on 
each  interval  and  each  basis  function  is  non-zero  on  no  more 
than  two  intervals.   Thus  on  each  interval  the  quadrature  is 
taken  with  the  functions  and  each  of  the  four  non-zero  basis 
functions. 

For  example,  on  the  I    x-interval  and  the  J    y-interval, 
we  evaluate  the  function  at  the  nine  quadrature  points  and 
take  its  product  with  the  value  of  the  basis  functions  at 
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the  point  and  the  appropriate  Gauss  weights,  and  then  sum 
over  the  nine  points  in  the  interval  for  each  of  the  four 
basis  functions.   This  procedure  is  continued  by  moving  to 
the  next  interval,  and  using  the  appropriate  basis  functions 
for  that  interval. 
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